#ifndef ARRAY_EVALUATOR
#define ARRAY_EVALUATOR

namespace array_helper
{

template <class operation>
struct array_evaluator
{
	template <class array_type, class expr_type>
	inline static void evaluate( array_type&& arr, expr_type&& expr )
	{
		typedef typename std::remove_reference<array_type>::type array_plain_type;

	    if ( arr.num_elements() == 0 ) return;

        if ( array_plain_type::n_dim == 1 )
            return evaluate_stack_traversal_1( arr, expr );
        else
            return evaluate_stack_traversal_N( arr, expr );
	}
	
	template <class array_type, class expr_type>
	inline static void evaluate_stack_traversal_1( array_type&& arr, expr_type&& expr )
	{
		typedef typename std::remove_reference<array_type>::type array_plain_type;
		static constexpr int first_dim = 0;
		
		auto iter = arr.begin();
		iter.load_stride(first_dim);
		expr.load_stride(first_dim);

		bool use_unit_stride = 
			iter.is_unit_stride(first_dim) && 
			expr.is_unit_stride(first_dim);

		int common_stride = std::max( expr.suggest_stride(first_dim), iter.suggest_stride(first_dim) );
		bool use_common_stride = 
			iter.is_stride( first_dim, common_stride ) && 
			expr.is_stride( first_dim, common_stride );

		typedef typename array_plain_type::num_type num_type;
		const num_type* last = iter.data() + arr.shape_[first_dim] * arr.stride_[first_dim];

		if( use_unit_stride || use_common_stride )
		{
			int ubound = arr.shape_[first_dim] * common_stride;
			num_type* data = iter.data();

			if( common_stride == 1 )
			{
				int n1 = ubound & 3;
				int i = 0;
				for(; i < n1; ++i )
				{
					const auto expr_data = expr.fast_read(i);
					operation::apply( *data++, expr_data );
				}
				
				for(; i < ubound; i += 4 )
				{
					const int t1 = i + 1;
					const int t2 = i + 2;
					const int t3 = i + 3;
					const auto tmp1 = expr.fast_read(i);
					const auto tmp2 = expr.fast_read(t1);
					const auto tmp3 = expr.fast_read(t2);
					const auto tmp4 = expr.fast_read(t3);
					operation::apply( *data++, tmp1 );
					operation::apply( *data++, tmp2 );
					operation::apply( *data++, tmp3 );
					operation::apply( *data++, tmp4 );
				}
			}
			else 
			{
				int n1 = ( arr.shape_[first_dim] & 3 ) * common_stride;

				int i = 0;
				for(; i != n1; i += common_stride )
				{
					const auto expr_data = expr.fast_read(i);
					operation::apply( data[i], expr_data );
				}
				
				int stride_inc = 4 * common_stride;
				for(; i != ubound; i += stride_inc )
				{
					const int i2 = i + common_stride;
					const int i3 = i + 2 * common_stride;
					const int i4 = i + 3 * common_stride;
					const auto tmp1 = expr.fast_read(i);
					const auto tmp2 = expr.fast_read(i2);
					const auto tmp3 = expr.fast_read(i3);
					const auto tmp4 = expr.fast_read(i4);					
					operation::apply( data[i],  tmp1 );
					operation::apply( data[i2], tmp2 );
					operation::apply( data[i3], tmp3 );
					operation::apply( data[i4], tmp4 );
				}
			}
		}
		else 
		{
			while( iter.data() != last )
			{
				operation::apply( *iter, *expr );
				iter.advance();
				expr.advance();
			}
		}

	}
	
	template <class array_type, class expr_type>
	inline static void evaluate_stack_traversal_N( array_type&& arr, expr_type&& expr )
	{
		typedef typename std::remove_reference<array_type>::type array_plain_type;
		
		/*
		 * A stack traversal replaces the usual nested loops:
		 *
		 * for (int i=A.lbound(firstDim); i <= A.ubound(firstDim); ++i)
		 *   for (int j=A.lbound(secondDim); j <= A.ubound(secondDim); ++j)
		 *     for (int k=A.lbound(thirdDim); k <= A.ubound(thirdDim); ++k)
		 *       A(i,j,k) = 0;
		 *
		 * with a stack data structure.  The stack allows this single
		 * routine to replace any number of nested loops.
		 *
		 * For each dimension (loop), these quantities are needed:
		 * - a pointer to the first element encountered in the loop
		 * - the stride associated with the dimension/loop
		 * - a pointer to the last element encountered in the loop
		 *
		 * The basic idea is that entering each loop is a "push" onto the
		 * stack, and exiting each loop is a "pop".  In practice, this
		 * routine treats accesses the stack in a random-access way,
		 * which confuses the picture a bit.  But conceptually, that's
		 * what is going on.
		 */

		/*
		 * ordering(0) gives the dimension associated with the smallest
		 * stride (usually; the exceptions have to do with subarrays and
		 * are uninteresting).  We call this dimension max_dim; it will
		 * become the innermost "loop".
		 *
		 * Ordering the loops from ordering(N_rank-1) down to
		 * ordering(0) ensures that the largest stride is associated
		 * with the outermost loop, and the smallest stride with the
		 * innermost.  This is critical for good performance on
		 * cached machines.
		 */
		const int max_dim = arr.storage_.get_storage_order(0);

		// Create an iterator for the array receiving the result
		auto iter = arr.begin();

		// Set the initial stack configuration by pushing the pointer
		// to the first element of the array onto the stack N times.
		int i;
		for ( i = 1; i < array_plain_type::n_dim; ++i )
		{
			iter.push(i);
			expr.push(i);
		}

		// Load the strides associated with the innermost loop.
		iter.load_stride(max_dim);
		expr.load_stride(max_dim);

		/* 
		 * Is the stride in the innermost loop equal to 1?  If so,
		 * we might take advantage of this and generate more
		 * efficient code.
		 */
		bool use_unit_stride = 
			iter.is_unit_stride(max_dim) && 
			expr.is_unit_stride(max_dim);
		
		/*
		 * Do all array operands share a common stride in the innermost
		 * loop?  If so, we can generate more efficient code (but only
		 * if this optimization has been enabled).
		 */
		int common_stride = std::max( expr.suggest_stride(max_dim), iter.suggest_stride(max_dim) );
		bool use_common_stride = 
			iter.is_stride( max_dim, common_stride ) &&
			expr.is_stride( max_dim, common_stride );

		/*
		 * The "last" array contains a pointer to the last element
		 * encountered in each "loop".
		 */
		typedef typename array_plain_type::num_type num_type;
		const num_type* last[ array_plain_type::n_dim ];

		// Set up the initial state of the "last" array
		for ( i = 1; i < array_plain_type::n_dim; ++i )
		{
			auto order_i = arr.storage_.get_storage_order(i);
		    last[i] = iter.data() + arr.shape_[order_i] * arr.stride_[order_i];
		}

		int last_length = arr.shape_[max_dim];
		int first_non_collapsed_loop = 1;

		/*
		 * This bit of code handles collapsing loops.  When possible,
		 * the N nested loops are converted into a single loop (basically,
		 * the N-dimensional array is treated as a long vector).
		 * This is important for cases where the length of the innermost
		 * loop is very small, for example a 100x100x3 array.
		 * If this code can't collapse all the loops into a single loop,
		 * it will collapse as many loops as possible starting from the
		 * innermost and working out.
		 */

		// Collapse loops when possible
		for ( i = 1; i < array_plain_type::n_dim; ++i)
		{
			// Figure out which pair of loops we are considering combining.
			int outer_loop_dim = arr.storage_.get_storage_order(i);
			int inner_loop_dim = arr.storage_.get_storage_order(i-1);

			/*
			 * The can_collapse() routines look at the strides and extents
			 * of the loops, and determine if they can be combined into
			 * one loop.
			 */

			if ( iter.can_collapse( outer_loop_dim,inner_loop_dim ) 
			  && expr.can_collapse( outer_loop_dim,inner_loop_dim ) )
			{
				last_length *= arr.shape_( outer_loop_dim );
				first_non_collapsed_loop = i + 1;
			}
			else  
				break;
		}
		
		/*
		 * Now we actually perform the loops.  This while loop contains
		 * two parts: first, the innermost loop is performed.  Then we
		 * exit the loop, and pop our way down the stack until we find
		 * a loop that isn't completed.  We then restart the inner loops
		 * and push them onto the stack.
		 */
		while (true) 
		{
			/*
			* This bit of code handles the innermost loop.  It would look
			* a lot simpler if it weren't for unit stride and common stride
			* optimizations; these clutter up the code with multiple versions.
			*/

			if ((use_unit_stride) || (use_common_stride))
			{
				/*
				 * The check for BZ_USE_FAST_READ_ARRAY_EXPR can probably
				 * be taken out.  This was put in place while the unit stride/
				 * common stride optimizations were being implemented and
				 * tested.
				 */

				// Calculate the end of the innermost loop
				int ubound = last_length * common_stride;

				/*
				 * This is a real kludge.  I didn't want to have to write
				 * a const and non-const version of FastArrayIterator, so I use a
				 * const iterator and cast away const.  This could
				 * probably be avoided with some trick, but the whole routine
				 * is ugly, so why bother.
				 */
				num_type* data = iter.data();

				/*
				 * BZ_NEEDS_WORK-- need to implement optional unrolling.
				 */
				if( common_stride == 1 )
				{
					/*
					for( int i = 0; i < ubound; ++i )
					{
						const auto& expr_data = expr.fast_read(i);
						operation::apply( data++, &expr_data );
					}
					*/
					auto task = [&data, &expr]( int i )
					{
						const auto expr_data = expr.fast_read(i);
						operation::apply( *data++, expr_data );
					};
					unroll_loop<4>::loop( ubound, task );
				}
				else 
				{
					/*
					for( int i = 0; i != ubound; i += common_stride )
					{
						const auto& expr_data = expr.fast_read(i);
						operation::apply( &(data[i]), &expr_data );
					}
					*/
					auto task = [&data, &expr, &common_stride]( int i )
					{
						const auto expr_data = expr.fast_read( i * common_stride );
						operation::apply( data[ i * common_stride ], expr_data );
					};
					unroll_loop<4>::loop( ubound / common_stride, task );
				}
				/*
				 * Tidy up for the fact that we haven't actually been
				 * incrementing the iterators in the innermost loop, by
				 * faking it afterward.
				 */
				iter.advance( last_length * common_stride );
				expr.advance( last_length * common_stride );
			}
			else 
			{
				/*
				 * We don't have a unit stride or common stride in the innermost
				 * loop.  This is going to hurt performance.  Luckily 95% of
				 * the time, we hit the cases above.
				 */
				typename array_plain_type::num_type* end = 
					iter.data()	+ last_length * arr.stride_(max_dim);

				while( iter.data() != end )
				{
					operation::apply( *iter, *expr );
					iter.advance();
					expr.advance();
				}
			}


			/*
			 * We just finished the innermost loop.  Now we pop our way down
			 * the stack, until we hit a loop that hasn't completed yet.
			 */ 
			int j = first_non_collapsed_loop;
			for(; j < array_plain_type::n_dim; ++j )
			{
				// Get the next loop
				int r = arr.storage_.get_storage_order(j);

				// Pop-- this restores the data pointers to the first element
				// encountered in the loop.
				iter.pop(j);
				expr.pop(j);

				// Load the stride associated with this loop, and increment
				// once.
				iter.load_stride(r);
				expr.load_stride(r);
				iter.advance();
				expr.advance();

			    // If we aren't at the end of this loop, then stop popping.
			    if( iter.data() != last[j] )
					break;
			}

			// Are we completely done?
			if( j == array_plain_type::n_dim )
				break;

			// No, so push all the inner loops back onto the stack.
			for( ; j >= first_non_collapsed_loop; --j )
			{
				int r2 = arr.storage_.get_storage_order(j-1);
				iter.push(j);
				expr.push(j);
				last[j-1] = iter.data() + arr.shape_[r2] * arr.stride_[r2];
			}

			// Load the stride for the innermost loop again.
			iter.load_stride( max_dim );
			expr.load_stride( max_dim );
		} //!< while(true)	
	}	
};

} //!< array_helper

#endif //!< ARRAY_EVALUATOR